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^T)". Abstract. The unified gas kinetic scheme (UGKS) of K. Xu et al. [M], originally 

developed for multiscale gas dynamics problems, is applied in this paper to a linear kinetic 
\ model of radiative transfer theory While such problems exhibit purely diffusive behavior 
in the optically thick (or small Knudsen) regime, we prove that UGKS is still asymptotic 
preserving (AP) in this regime, but for the free transport regime as well. Moreover, this 
scheme is modified to include a time implicit discretization of the limit diffusion equation, 
and to correctly capture the solution in case of boundary layers. Contrary to many AP 
schemes, this method is based on a standard finite volume approach, it does neither use any 
<^ ' decomposition of the solution, nor staggered grids. Several numerical tests demonstrate the 
properties of the scheme. 
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1 Introduction 

Kinetic models are efficient tools to describe the dynamics of systems of particles, like in 
rarefied gas dynamics (RGD), neutron transport, semi-conductors, or radiative transfer. 
Numerical simulations based on these models require important computational resources, 
q ! but modern computers make it possible to simulate realistic problems. 

These simulations can be made much faster when the ratio between the mean free path of 
particles and a characteristic macroscopic length (the so-called Knudsen number in RGD, de- 
noted by e in this paper) is small. In such cases, the system of particles is accurately described 
by a macroscopic model (Euler or Navier-Stokes equations in RGD, diffusion equations in 
neutron or photon transport) that can be numerically solved much faster than with kinetic 
models. 

However, there are still important problems in which the numerical simulation is difficult: 
in multiscale problems, e can be very small in some zones, and very large elsewhere (opaque 
vs. transparent regions in radiative transfer). Standard numerical methods for kinetic equa- 
tions are very expensive in such cases, since, for stability and accuracy reasons, they must 
resolve the smallest microscopic scale, which is computationally expensive in small e zones. 
By contrast, macroscopic solvers are faster but may be inaccurate in large e zones. 
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This is why multiscale numerical methods have been presented in the past 20 years: the 
asymptotic-preserving (AP) schemes. Such schemes are uniformly stable with respect to e 
(thus their computational complexity does not depend on e), and are consistent with the 
macroscopic model when e goes to (the limit of the scheme is a scheme for the macroscopic 
model). 

AP schemes have first been studied (for steady problems) in neutron transport by Larsen, 
Morel and Miller [25], Larsen and Morel [24] , and then by Jin and Levermore [T2| IT5] . For 
unstationary problems, the difficulty is the time stiffness due to the collision operator. To 
avoid the use of expensive fully implicit schemes, several semi-implicit time discretizations 
schemes, based on a decomposition of the distribution function between an equilibrium part 
and its deviation, have been proposed by Klar [19], and Jin, Pareschi and Toscani [18J (see 
preliminary works in [T7| [TTj and extensions in JT6J [151 EB [20j [21]). Similar ideas have also 
been used by Buet et al. in [2], Klar and Schmeiser [22], Lemou and Mieussens [281 H]> 
and Carrillo et al. [U [5] . The theory of well balanced schemes is another way to obtain AP 
schemes, as in the work of Gosse an Toscani JS[ [S]. Other approaches have been recently 
proposed by Lafitte and Samaey [23] and Gosse [7], but there extensions to more complex 
cases is not clear. Finally, the idea of [H] has been renewed to obtain an AP scheme for 
linear equations on two-dimensional unstructured meshes in the work of Buet, Despres, and 
Franck [5J. All these methods have advantages and drawbacks, and there is still a need for 
other AP schemes. 

A rather different approach has recently been proposed by K. Xu and his collaborators, 
in the context of rarefied gas dynamics [54] . This method is called unified gas kinetic 
scheme (UGKS) and is based on a gas kinetic scheme which has been developed by K. 
Xu since 2000 (see |33j for the first reference and many other references in [M]). Roughly 
speaking, the UGKS is based on a finite volume approach in which the numerical fluxes 
contain information from the collision operator. In some sense, it has some connexions with 
the well balanced schemes developed for hyperbolic problems with source terms in [14] [HI [9] , 
even if the construction is completely different. While this approach to design AP schemes 
looks very promising, it seems that it has not yet received the attention it deserves from the 
kinetic community. This is probably due to the fact that the nice properties of the UGKS 
presented in [5J] are difficult to understand for people who are not specialist of gas kinetic 
schemes. 

However, we believe that the UGKS approach is very general and can benefit to many 
different kinetic problems. Let us mention that the big advantage of the UGKS with respect 
to other methods is that is does not require any decomposition of the distribution function 
(hence there is no problem of approximation of the boundary conditions), it does not use 
staggered grids (which is simpler for multi-dimensional problems), and it is a finite volume 
method (there is no need of discontinuous Galerkin schemes that are more expensive). 

In this paper, our first goal is to present UGKS in a very simple framework, so that it 
can be understood by any researcher interested in numerical method for kinetic equations. 
We also want to show that the UGKS can be successfully applied to other fields than RGD. 
Here, it is used to design an AP scheme for linear kinetic equations, namely a simple model 
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of radiative transfer. Such an extension is not obvious, since linear models exhibit a purely 
diffusive (parabolic) behavior in the small e regimes, while models from RGD (like the 
Boltzmann equations) have a rather convection (hyperbolic) behavior. Indeed, even if the 
UGKS is originally made to correctly describe this convection regime and to capture the small 
viscous effects (like in the compressible Navier-Stokes equations), we prove in this paper that 
it can also capture a purely diffusive effect. Moreover, we propose several extensions: implicit 
diffusion, correct boundary conditions for boundary layers, treatment of collision operator 
with non isotropic scattering kernel. The scheme is proved to be AP in both free transport 
and diffusion regimes, and is validated with several numerical tests. 

The outline of our paper is the following. In section EJ we present the linear kinetic model, 
and its approximation by the UGKS. Its asymptotic properties are analyzed in section [3j 
Some extensions are given in section HJ arid the scheme is validated with various numerical 
tests in section [5j 



2 The UGKS for a linear transport model 

2.1 A linear transport model and its diffusion limit 

The linear transport equation is a model for the evolution (by transport and interaction) of 
particles in some medium. In this paper, we are mainly concerned by the radiative transfer 
equation, which reads, in scaled variables 

ed t <b + Q ■ V r (b = -(— [ cbdQ-d)) - each + eG, 
e 4tt J 

where (p(t, r, fl) is the spectral intensity in the position-direction phase space that depends 
on time t, position r = (x, y, z) e M 3 , and angular direction of propagation of particles 
Q G S 2 . Moreover, a is the scattering cross section, a is the absorption cross section, and 
G is an internal source of particles. These three last quantities may depend on x, but they 
are independent of Q. The linear operator i-> j- J (pdfl — models the scattering of the 
particles by the medium and acts only on the angular dependence of <fi. This simple model 
does not allow for particles of possibly different energy (or frequency); it is called "one- 
group" or "monoenergetic" equation. The parameter e is a scale factor that measures the 
ratio between a typical microscopic length (the mean free path of a particle, for instance) to 
a typical macroscopic length (the size of the computational domain, for instance). See [25] 
for details. 

In this paper, we consider this one-group equation in the slab geometry: we assume that 
4> depends only on the slab axis variable x 6 R. Then it can be shown that the average 
of (j) with respect to the (y,z) cosine directions of Q, denoted by f(t,x,v), satisfies the 
one-dimensional equation 

ed t f + vdj = -((f) - f) - eaf + eG, (1) 
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where v G [—1,1] is the x cosine direction of Q and the operator (.) is such that ((f)) = 
\ J"* (p{v) dv is the average of every f-dependent function <p. 

When e becomes small, it is well known that the solution / of (CQ) tends to its own average 
density p = (f), which is a solution of the asymptotic diffusion limit 

d t p + d x nd x p = -ap + G, (2) 

where the diffusion coefficient is k(x) = = ~~ 3(T ^ ■ An asymptotic preserving scheme 

for the linear kinetic equation ([I]) is a numerical scheme that discretizes ([I]) in such a way 
that it leads to a correct discretization of the diffusion limit (jSJ) when e is small. 



2.2 First ingredient of the UGKS: a finite volume scheme 

The first ingredient of the UGKS is a finite volume approach. Equation ([1]) is integrated over 
a time interval [t n ,t n+ i] and over a space cell [x i+ i,Xj_i] to obtain the following relation: 

f n+1 - f n 1 / \ 1 f tn+1 f x i+i a 

where /" = ^ f{t n ,x,v)dv is the average of / over a space cell and is the 
microscopic flux across the interface ^.i: 
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eAt 



vf(t,x i+ i,v)dt. (3) 



To obtain a scheme which is uniformly stable with respect to e, the collision term, which is 
the stiffest term in the previous relation when e is small, must be discretized by an implicit 
approximation. For simplicity, we use here the standard right-rectangle quadrature which 
gives a first order in time approximation: 

i rtn+i r x i+ i rt 



AtAx 



/ / ^-(p-n-afjdxdt^^ipr'-fr 1 )-^ 

Jtn Jx. i E 6 

1 2 



We also assume that the scattering and absorption coefficient do not vary to much inside 
a cell, so that the average of the products af and af are close to the product of the 
averages of each terms. Taking the absorption term implicit is not necessary, and an explicit 
approximation could also be used. 

Then the finite volume scheme for (prj) reads: 



f n+l _ fn 1 . \ (J 

i% jt | ( i i \ _ u% („n+l fn+l\ „, fn+l 



+ tt <t> i+ ± - = 3(pr +i - fri - <xifr l + g. (4) 



Now, it remains to approximate the flux , and hence the value of / at any time t between 
t n and t n+ i at the cell interface x i+ i, by using the averaged values f™ +1 , f™_ X) etc. The 
choice made for this approximation is the core of the UGKS scheme. 
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Note that in this construction, v is not discretized: indeed, while it can be discretized 
by any method, this would not change our analysis. Keeping v continuous at this stage 
considerably simplifies the notations, what we do here. See section |5] for an example of 
velocity discretization. 



2.3 Second ingredient: a characteristic based approximation of 
the cell interface value f(t,x i+ i,v) 



In this section, we explain how f(t,x i+ i,v) is reconstructed in the flux defined in To 
emphasize the importance of this approximation, note that the standard first-order upwind 
approximation 

f{t,x i+ i,v) = f?l v>0 + /r +1 i,<o, 

where l^o = 1 if v 5s and else, is not a choice that gives an AP scheme. Indeed, it 
introduces numerical dissipation that dominates the physical diffusion in the diffusive limit 
regime. This can also be interpreted as follows: this approximation is nothing but the 
solution at time t of the Riemann problem 

d t f + -vdj = 0, 

£ 

f(tn, x, v) = f" if x > x i+ i , f? +1 else. 

This problem does not take into account the right-hand side of the real problem ([1]). It is 
well known that hyperbolic problem with source term cannot be accurately approximated if 
the numerical flux is constructed by ignoring the source term, in particular in limit regimes 
(see [H], for instance). While it is not always easy to take the source term into account in 
the numerical flux, the transport+relaxation structure of ([I]) makes this problem particularly 
simple. 

Indeed, in [34], Xu and Huang propose to use the integral representation of the solution 
of the BGK equation, which, for our model (JTJ, is obtained as follows. If the coefficients a 
and a are constant in space and time, equation (CQ) is equivalent to 

-e"7(t, x + -t, v) = e ut -pit, x + -t)+i 
at e \e z e 

where v — ^ + a. If they are not constant but slowly varying in one cell, we consider this 

relation as an approximation around t n and x i+ i, and we denote by cr i+ i, a i+ i and u i+ i the 

corresponding constant values of a, a and v. Then we integrate this relation between t n and 

some t < t n+ i, and we replace x + |£ by x i+ i to get the following relation: 

f(t,x i+l v) « e' V ^ {t ' tn) f{t n ,x i+k - -(t-t n )) 

2 2 £ 

+ I e~ U ^ (t - tn) ^ P (s, x i+ , - V -(t - s)) ds (5) 
l_ e -"« + j ('-*"> 
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Now, it remains to design an approximation of two terms: the first one is / at t n around 
that is to say f(t n ,x i+ i — ^(t — t n )), and the second one is p between t n and t n+ i 
around x i+ i, that is to say p(s,x i+ i — ^(t — s)). 

To approximate / at t n around x i+ i, the simplest approach is to use a piecewise constant 
reconstruction: 

f™ if x < Xa , 1 
f(t n ,x,v) = { +5 (6) 

if x > X l+ i 
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Of course, a more accurate reconstruction can be obtained (for instance a piecewise linear 
reconstruction with slope limiters, as in [31]), but the presentation of the scheme is simpler 
with this zeroth order reconstruction. 

This is for the approximation of p between t n and t n+ \ around x i+ i that we need the 
second main idea of K. Xu: this reconstruction is piecewise continuous. This can be sur- 
prising, since p is the velocity average of / which is represented by piecewise discontinuous 
function, but this is the key idea that allows the scheme to capture the correct diffusion 
terms in the small e limit. In [34j HO] , this reconstruction is piecewise linear in space and 
time. However, in the context of the diffusion limit, we found that the a piecewise constant in 
time reconstruction is sufficient, while a piecewise continuous linear reconstruction in space 
is necessary. First, we define the unique interface value p n , at time t n by: 



Pr + i = </fl,>o + /f + il,<o>- (7) 
Then, we define the following reconstruction for t in [t n ,t n+ i] and x around x i+ i: 



P™+i + $xp";+i (x - x i+ i) if x < x i+ i 



p{t,x) 



P n i+ l + $xP" + i(x ~ x i+ i) ifx>x i+ i 
with left and right one-sided finite differences: 

P i+ L Pi D Pi+1 Fj+i 

Remark 2.1. When e is very small, the foot of the characteristic x i+ i — ~(i — s) in (JSJ) 
can be very far from x i+ i, and hence using the reconstruction (jHJ) might be very inaccurate. 
However, note that in flSj), p is multiplied by an exponential term which is very small in 
this case, and we can hope that the inaccuracy made in the reconstruction (jSJ) has not a 
too strong influence. Indeed, our asymptotic analysis and our numerical tests show that the 
accuracy of the scheme is excellent, see sections [3] and |5j 
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2.4 Numerical flux 

Now, the numerical flux (j> i+ i = ^ v f{^i x i+i, v ) dt can be computed exactly by using 
expressions ( I5l6ll8j) to get 



i =A+^(/ri,>o+/r+ii,<o) 



(10) 
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■J<0j 



where the coefficients = A(At, e, cr i+ i , a i+ i), C i+ i = C(At, e, cr i+ i , a i+ i ), A+i = 

D(At, e, CKj+i), and -Ej+i = E(At, e, a i+ i,a i+ i) are defined by the following functions: 

1 



A(At,e,a,a) = ^ ; (l-e-^) (11) 



C(At, e, *,a) = -^ (At - i (l - e^*) ) (12) 
D(At,.,a,a) = -^(Af (1 + e^) - \(l - e"^)) (13) 
25(At, e , ff) a ) = ^_ (A* - 1 (1 - e^) ) (14) 



where we remind that z/ = + a. 

Now, using (j3J), /™ + can be o 
is done by using the conservation law in the following section. 



£ 

Now, using P|, can be obtained, provided that we can first determine This 



2.5 Conservation law 

It is now well known that semi-implicit schemes for relaxation kinetic equations can be solved 
explicitly by using the corresponding discrete conservation laws: see [31] where it was first 
suggested for the BGK equation, and [34] for the use of this technique to design the UGKS. 

The idea is to eliminate in (j3J) by taking its f-average, hence obtaining the following 
discrete conservation law 

' ; % + h ~ = -aiP? +1 + G (15) 



At ' Ax V~ i+ i ~*-f 
where the macroscopic numerical flux is 

By using (flOl . we find: 



2 



1 n n — rf< 



Note that relation (fl~5l) is still implicit, but it can be solved explicitly. 
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2.6 Summary of the numerical scheme 

Finally, is computed as follows for every cell i: 
1. compute by solving: 



At Ax 
2. compute by solving: 



At As V l+ 2 e 2 

where macroscopic and microscopic fluxes are given by (|T6j) and (TlOT) . 

Of course, this scheme must be supplemented by numerical boundary conditions. This 
will be detailed in the next section, after the asymptotic analysis. 

3 Asymptotic analysis 
3.1 Free transport regime 

The behavior of the scheme in the small a, a limit is completely determined by the following 
property of the coefficient functions A, C, D, and E. 

Proposition 3.1. When a and a tend to (while e is constant), we have: 

• A(At, e, a, a) tends to - 

• C and D(At,e,a,a) tend to 

• E(At,e,a,a) tends to 

As a consequence, the microscopic flux <f) i+ ± defined in ( ITUl) has the following limit: 

—7> ~07%>o + /r +1 l,<o) + ^vG. 

2 a— i>U S L 

This is nothing but a consistent first-order upwind flux (plus a constant term that has no 
influence in the scheme), and the limit finite volume scheme (j4]) is 

rn+l _ rn i 

At + Ai~e ((//%>0 + f ^ +llv<o) ~ (/ -^>° + ^^< o) ) = G ' 

which is indeed a consistent approximation of the limit transport equation ([T|) when a and 
a tend to 0. As a consequence, the UGKS is AP in this limit. 
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3.2 Diffusion regime 



Similarly, the behavior of the scheme in the small e limit is completely determined by another 
property of the coefficient functions A and D. 

Proposition 3.2. When e tends to 0, we have: 

• A(At, e, a, a) tends to 

• D(At,e,a,a) tends to —-. 

As a consequence, the macroscopic flux defined in f|T6|) has the following limit: 

— -* - Q Pi+ \ ~ , (17) 

i+ 2 

and hence the limit of the discrete conservation law (1151) is 



Pi Pj_ 1 / 1 Pi+1 Pi 1 Pi 

At Ax V 3cr,- , i Ax 3ct,;_ i Ax 
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which is a consistent approximation of the diffusion limit (J2J), with a standard three-point 
centered approximation of the second order derivative of p. This proves that UGKS is AP 
for the limit e — > 0. 



3.3 Boundary conditions 

If we consider equation ([T]) for x in the bounded domain [0, 1], we need the following boundary 
conditions: 

f(t,x = 0,v>0) = f L (v) and f(t,x = l,v<0) = f R (v), (19) 

where fi and fn can be used to model inflow or reflexion boundary conditions. In the 
case of an inflow boundary condition, if fi and f R are independent of v, the diffusion limit 
is still (j2J) with the corresponding Dirichlet boundary data p(t, 0) = fz and p(t, 1) = f R . 
However, when one of the boundary data (say fi) depends on v (which is called a non 
isotropic data), using this diffusion approximation requires some modifications. First, a 
boundary layer corrector must be added to p to correctly approximate /. Moreover, / is well 
approximated by the solution p of the diffusion equation outside the boundary layer only if 
this equation is supplemented by the boundary condition 

p(t, x = 0) = / W(v)f L (v) dv = 2(Wf L l v>0 ) (20) 

where W(v) is a special function that can be well approximated by 0.956v + 1.565f 2 ~ 
3/2v 2 + v (see [25]). 
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Now, we explain how these boundary conditions can be taken into account numerically in 
the UGKS. We assume that we have i max space cells in [0, 1]. As compared to our derivation 
of scheme (j4]), the difference is that the derivation of the numerical fluxes at the boundaries 
07 and (f) n i must take into account the boundary data. For simplicity, we only describe 

in details how the left boundary flux is constructed. The definition §5§ is still valid, but the 
integral representation of / at this left boundary now is 

r h if ^ > o 

e tn) f(t n ,xi - -(t - t n )) + / e *V tn) ^p(s,xi--(t-s))ds 

— Vl(t—tn) 

1 — e 2 

+ G ifv<0 

2 



f(t,xi,v) 



where f(t n ,xi — ^(t — t n )) and p(s,xi — -(t — s)) have to be defined for v < only, that 
is to say on the right-hand side of the left boundary. According to the approximations ([6]) 
and (jSJ), we set for v < and x > xi: 

f(t n ,x,v) = 

and for t in [t n ,t n+ i] 



P{t,x) =p n i+S x pT(x-xi), 



2 

where according to (JTJ) the left boundary value of p now should be 



p^H/ii^o + zrw, (21) 

2 

Pl-P\ 

and 5 x p\ = Aa , . With this definitions, the numerical flux at the left boundary is found 
to be 

0i = -hl v >o + ^i^/i n l,<o + Civpl l v<0 + DivXp^l^o + EivGl v<0 , (22) 

25 2 s - 25 22 2 

where the coefficients Ai, Bi, Ci, Di and Ei are defined as in section l2~4l The corresponding 

2 2 2 22 

macroscopic flux is 

$1 = -{vf L lv>o)+Ai (vf?t v<0 ) + (vl v<0 )Cxpl + (v 2 l v<0 )Di5 x pl R l v<0 + (vl v<0 )EiGl v<0 . 

(23) 

We could have used the exact values (vt v<0 ) = —1/4 and (f 2 l„ <0 ) = 1/6 in ([23]) . but in 

practice, when v is discretized, it is better to use the numerical approximation of (vt v<0 ) 

and (v 2 t v< o) corresponding to the chosen quadrature formula: this avoids an important loss 

of accuracy when e goes to 0. 

With these relations that define 61 and $1, the numerical scheme summarized in sec- 
2 2 

tion 12.61 can be used in every cell i from 1 to i max . 
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Unfortunately, this scheme is not uniformly stable for small e, hence cannot be AP in 
the diffusion limit. The reason is that the macroscopic flux now contains two unbounded 
terms that are ^(vfil v> o) and — \C\_p\ (we can prove that C is asymptotically equivalent 

to -). This fact is not observed in [34], since UGKS is not applied in a diffusive scaling, and 
hence the boundary data does not contribute as a 1/e term in the boundary flux. However, 
this drawback can be easily corrected. First, we propose the following simple correction: it 
is sufficient to modify the definition (1211 of the boundary value p\ to 

2 

A = -4^, (24) 

so that the two unbounded terms exactly cancel each other. This definition ensures that 
$i is uniformly bounded with respect to e, and it is reasonable to conjecture that the 
corresponding scheme is uniformly stable. 

Now, we investigate the two numerical limits of our scheme with this modified boundary 
condition. Note that this definition, while not consistent with the physical boundary data, 
is a standard approximation of the exact Dirichlet boundary condition of the diffusion limit. 
This means that we can hope for a correct behavior of the scheme in the diffusion limit. 
Indeed, following the same analysis as in section 13.21 we find that the discrete conservation 
law (|15|) in the first cell i — 1 tends to 

N 2 2 ' 

which is consistent approximation of the diffusion limit equation, with a Dirichlet boundary 
data p\ given by ( |24|) . This boundary data is the correct one if Jl is isotropic (since we 

2 

get p\ = Jl), but is only an (standard) approximation of the exact value ( |20l) if fa is not 

2 

isotropic. Note that a similar problem occurs with several AP schemes in the literature, like 
in the schemes of [T71 [28] , while it seems to be greatly reduced in the recent scheme of [27] . 
We propose in section I4T21 a modification of the boundary condition of our scheme to increase 
the accuracy near the boundaries. 

However, contrary to the schemes that have been mentioned, the inaccuracy of our mod- 
ified boundary value p\ has no influence in the free transport regime, which is another 

2 

interesting property of the UGKS. Indeed, since coefficients C, D, and E tend to for small 
a, the boundary fluxes <bi and $i tend to the simple first-order upwind boundary fluxes 

2 2 

4>i = -(f L lv>0 + fitv<o), 
2 E 

$1 = (0I> = -(Vfrtv>0 +Vf?t v<0 ). 

2 2 £ 

4 Extensions 

In this section, we propose some extensions of the previous scheme. 
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4.1 Implicit diffusion 

Our scheme gives for small e an explicit scheme for the limit diffusion equation. This 
means that at this limit, our scheme requires the following CFL like condition to be stable: 
At < where p = 1/ min(3cr) is the largest diffusion coefficient in the domain. When Ax 
is small, this condition is very restrictive. Indeed, diffusion equations are generally solved by 
implicit schemes that are free of such a CFL condition. It is therefore interesting to modify 
our scheme so as to recover an implicit scheme in the diffusion limit. Such a modification is 
not always trivial (it is not known for schemes of (TTJ |28j EIHJ El E3]), but has already been 
proposed for some others (see [221 126] ). 

For our scheme, we first note that the explicit diffusion term in ( [18]) comes from only one 
term in the microscopic numerical flux, namely the fourth term in (flQj) : D i+ iv 2 (5 x p n ^it v> o + 

5 x p n ^i l«<o)- Indeed, the left and right slopes 8 x p n ^ 1 and Sxp 71 ^ (defined by ([9])) depend on 

1+ 2 1+ 2 »+ 2 

the interface value p n , and on the left and right values p" and pf +1 at time t n . After 

l+ 2 

integration with respect to v, the interface value vanishes and we only get the difference 
of the pf and pf +1 at time t n . Consequently, to obtain an implicit diffusion term, that is 
to say a difference of left and right values of p at time t n +i, it is sufficient to modify the 
definition of the slopes as follows: we now set 

n _ „n+l n+1 _ n 

r Pi+l Pi D Pi+l 

«"tfi = ^fc7T- and «-^ = -S^' (26) 

Note that the interface value p n L is still at t n , so that the flux is still explicit with respect 

J+ 2 

to /. All the other quantities that are defined our scheme are unchanged. 

Now, following the same analysis as in section I3.2[ our scheme converges to the following 
implicit discrete diffusion equation: 

p- +1 -p? i ( i Ptl - pT x i pT x - ptl i „+, r 

— CtiPi "+■ Or. 



At Ax V 3a,,! Ax 3a, i Ax 



i+2 



4.2 More accurate boundary conditions for the diffusion limits 

Here, we propose a modification of the UGKS to obtain the correct Dirichlet boundary 
value (120]) in case of a non isotropic boundary condition. Since the incorrect value (|24p is 
required by the first term -{vfz,l v >o) of $i in (123]) . we propose to replace this term by the 
value that gives the correct limit. Indeed, we now impose the relation 

^ = -^^>(W/xW + ..., (27) 
where ■ • • stands for the other terms of (T2"3l that are unchanged. Moreover, we set 

pl = 2(Wf L l v>0 ) (28) 
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in (127)1 and (122|) . These definitions ensure that $1 is uniformly bounded with respect to e, 
and moreover, the discrete conservation law f|T5|) in the first cell i = 1 tends to (1251) with the 
correct Dirichlet boundary data 2(W /l1„>o)- 

However, contrary to the definition of p\ suggested in section 13.3) this definition has 

2 

an influence in the free transport regime. Indeed, while <pi tends to the correct first order 
upwind boundary flux, the macroscopic boundary flux $i tends to ^(W JlIvX) +vfit v< o), 
which is not the correct flux for the free transport regime (it should be - (v /l1 1 ,>o+ , w /flixo)) • 
Finally, we can obtain the two correct limits by using a blended modification of the first 
term of $i in fl23|) . Indeed, we set 

$i = -([(l-9(ui))v + 9(vi_) x (-2{vl v<0 ))W]f L l v>0 ) + --- , (29) 

where again ■ • • stands for the other terms of fl23l) . We set 

p| = 2( [((1 - 0(z/|))v + 0(vi)W] f L l v>0 ) (30) 

in ( 1271) and ( 1221) . and the blending parameter is #(^) = 1— exp(— uAt). This parameter is such 
that 9{v) tends to 1 for small e and tends to for small a and a. This definition implies that 
in the diffusion limit, we get indeed the correct Dirichlet boundary condition 2(W /l1«>o) , 
and that in the free transport regime, the microscopic and macroscopic boundary fluxes tend 
to the correct values. Note that these properties are illustrated in section [5j 

4.3 General linear collision operators 

In this section, we propose a way to adapt the UGKS approach to a general linear Boltzmann 
operator. Namely, we consider the equation 

ed t f + vd x f = hf, (31) 

where the collision operator now is 

Lf = J\v,v')(f(v')-f(v))dv'. (32) 

We ignore the absorption and source terms here, since this does not change our analysis. 
For such a model, it is well known that the density of / satisfies in the small e limit the 
following diffusion equation 

d t p + d x nd x p = 0, (33) 

with n = (t>L -1 t>), and where L~ l is the pseudo-inverse of L. This operator is defined for 
functions with zero average by the following: for any <p such that (0) = 0, ip — L~ l <p is the 
unique solution of Lip = (f> such that (ip) = 0. 
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Since in the UGKS approach the relaxation form of the collision operator is strongly 
used, a natural idea is to write Lf in the gain-loss form: Lf = L + f — ■^hxf, where 1/t(v) = 

f]_ x k{y,v') dv' . Then we can try to apply the previous strategy in which the gain term 
L + f plays the role of p. However, we observed that this strategy fails completely, since the 
resulting scheme cannot capture the correct diffusion limit. 

Then we propose a different and simple way to capture the correct asymptotic limit by 
using the penalization technique of Filbet and Jin [Bj. First, the collision operator is written 
as 

Lf = (Lf - ORf) + 9Rf, (34) 

where Rf — p — f (with p = (/)) is the corresponding isotropic operator, and where 9 is a 
parameter adjusted to capture the correct diffusion coefficient. Therefore, equation (l3~Ii) can 
be rewritten as 

ed t f + vd x f= 6 -(p-f)+eG, (35) 

where G = (Lf — 6Rf)/e 2 is considered as a time and velocity dependent source. Then our 
previous derivation can be readily applied to this equation: we replace G by G, o by 9, and 
a by in the different steps of sections 12.21 to 12.61 to get the following scheme: 

7 + 1 ) = 0, 

At Ax V J+ 5 

nn+1 _ rn 1 . \ 1 

The numerical fluxes are given by relations (|T0T[T4|) . in which we just have to replace G by 
G, a by 8, and a by (note that v must be modified accordingly in these relations). 

Consequently, our analysis detailed in section [3] directly applies to this scheme, and we 
obtain that it gives for small e this limit diffusion scheme 

P? +1 -P? 1 ( (v 2 )P? + i-P? , (v^ Pl-PU ^ _ Q 



At Ax V 9 Ax 9 Ax 

which is a consistent approximation of the diffusion limit (I3"3"|) . if and only if the parameter 
9 is set to 

»— 7#r- (36) 

This proves that this "penalized" UGKS is AP for the limit e — > 0. 

Note that to use this scheme, we just have to compute L _1 v, which has to be done 
only once. This computation can be made analytically (for simple scattering kernels) or 
numerically. Also note that an approach with similar aspects has been proposed in [26] . 

However, this AP property is true only if the scheme is uniformly stable. Since this 
property is difficult to obtain for the complete scheme, we restrict to the space homogeneous 
problem, and we show that this leads to a non trivial restriction on the collision operator. 
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Proposition 4.1. The penalized scheme 

-x P n+l -r +1 ) + -ALf n -6Rr) 



fn+l in n i 



that approximates the homogeneous equation d t f = Lf is absolutely stable if At^k^ — 9) < e 2 , 
where ku — max B y k. This stability is uniform with respect to e if 9 > kw 

This proposition can be easily proved by using the fact that if the scattering kernel is 
bounded (0 < k m < k(v,v') < k M ), then L is a non positive self-adjoint operator, and that 
its eigenvalues are all bounded in absolute value by 2/c M . 

Note that with 9 defined by ( 136]) . this restriction reads —(vL v) < (v 2 )/kM- If this 
inequality is not satisfied, then At must be smaller and smaller as e decreases, and hence 
the scheme cannot be AP. This property will be studied in detail in a forthcoming work for 
realistic functions k (like the Henyey-Greenstein function). 



5 Numerical results 

Almost all the test cases presented here are taken from references [P9] IP?] . Comparisons 
of several existing AP schemes with these test cases can be found in [28J. Depending on 
the regime, we compare UGKS to a standard upwind explicit discretization of ([1]) or to the 
explicit discretization of the diffusion limit fl2]). These reference solutions are obtained after 
a mesh convergence study: the number of points is sufficiently large to consider that the 
scheme has converged to the exact solution. Generally, the time step for the UGKS is taken 
as At = cfl max(eAx, 3Ax 2 a/2), where cfl = 0.9. 

5.1 Various test cases with isotropic boundary conditions 
Example 1 Kinetic regime: 

xe[o,i], h(v) = o, f R (v) = i, 

a = l, a = 0, G = 0, e = 1. 

The results are plotted at times t = 0.1, 0.4, 1.0, 1.6, and 4. We use 25 and 200 points for 
UGKS. The reference solution is obtained with 1000 points. In figure HJ we observe that the 
UGKS is very close to the reference solution, even with the coarse discretization for short 
times t = 0.1 and t = 0.4 (which is better than the AP schemes compared in [2"%]). 

Example 2 Diffusion regime: 

are [0,1], h(v) = l, f R {v) = 0, 
a = l, a = 0, G = 0, e=lQ- 8 . 
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The results are plotted at times t = 0.01, 0.05, 0.15 and 2. We use 25 and 200 points for 
the UGKS. Here, the reference solution is obtained with the explicit discretization of the 
diffusion equation, since the kinetic equation cannot be solved by a the standard upwind 
scheme with such a small e. In figure EJ we see that the UGKS and the diffusion solution 
are almost indistinguishable at any times for both coarse and fine discretizations. 

Example 3 Intermediate regime with a variable scattering frequency and a source term: 

xe[o,i], h(v) = o, f R (v) = o, 

a = 1 + (10x) 2 , a = 0, G = l, e = 10- 2 . 

The results are plotted at times t = 0.4 with 40 and 200 points in figure [3] The reference 
solution is obtained with the explicit scheme using 20 000 points. We observe that the UGKS 
provides results that are very close to the reference solution, like the AP schemes compared 
in [28]. 

5.2 Comparison of different numerical boundary conditions 

In this test, we compare the different numerical boundary conditions (BC) proposed for the 
UGKS in section [3731 and [4721 The BC given in ff24|) that ensures the stability of the scheme 
is called stabilized BC. We call corrected BC the one described in ( 12711281) that gives a correct 
boundary value in the diffusion limit. Finally, the blended BC is that defined in ( 12911301) . 

Example 4 First, we consider an intermediate regime with a non-isotropic boundary con- 
dition that generates a boundary layer of size L = 0.01 at the left boundary: 

XG[0,1], f L (v)=v, f R (v) = 0, 
a = l, a = 0, G = 0, e = 10~ 2 

The results are plotted at times t = 0.4 with 25 and 200 points in figure H] (see [28] for a 
comparison of other AP schemes on the same test case). The reference solution is obtained 
with the explicit scheme using 20 000 points. We also consider the diffusion scheme with 
25 points, with a left boundary condition pi = 17/24 computed by formula (|20|) with the 
approximation W(v) = 3/2v 2 + v. With the coarse discretization, the boundary layer is of 
course not resolved. The stabilized BC is rather different from the reference solution in all 
the domain. At the contrary, we observe that the corrected and blended BC are very close 
to the reference solution inside the domain (far from the boundary). 

With the fine discretization, the boundary layer is resolved: surprisingly, the corrected 
and blended BC are less close to the reference solution, while the stabilized version is now 
closer, even in the boundary layer. This fact is rather different from the other AP schemes 
already used for the this test (see [28j). We do not understand well the reason so far, but it 
might be linked to the fact that the coefficients A, C, D of the UGKS depend in an intricate 
way on the numerical parameter At, Ax and e: for instance for a given small e, A and D 
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are closer to their limit value and —1/a for large At and Ax than for smaller values. See a 
similar observation in [32] about the high accuracy of the gas kinetic scheme (a macroscopic 
version of the UGKS) for under resolved cases. 

Example 5 Then we take e = 10~ 4 to be in a fully diffusion regime. The boundary layer 
is here very small and cannot be captured. The reference solution then is the solution of the 
diffusion equation, computed with 2000 cells. The results are plotted at t — 0.4 in figure El 
The comparison is almost the same as with e = 10~ 2 for the coarse mesh: the stabilized BC 
gives a result which is not very accurate, while, as expected, the corrected and blended BC 
(which give the same results) are much closer to the diffusion solution. However, with the 
fine mesh, the corrected and blended BC now give no difference with the diffusion solution, 
while the stabilized BC is not accurate. 

Example 6 Finally, in order to test our blended BC in all the regimes, we take e — 1 and 
a = 1, which is a purely transport regime (no collision). The results are plotted at times 
t = 0.4 with 25 and 200 points in figure [6j Here, the comparison with a reference solution is 
difficult: indeed, it is well known that discrete ordinates methods in free transport regimes 
require many velocity points to be accurate. With a small number of velocities (like the 
16 Gauss points we use here), the results show the well known "ray effect" that looks like 
several plateaus in the solution, and this phenomenon is stronger when the space mesh is 
refined. Consequently, we think it more interesting to compare the UGKS to a standard 
upwind explicit scheme with the same space resolution. As expected, we observe in figure [6] 
that both stabilized and blended BCs are very close to the standard scheme (for both coarse 
and refine space grids). At the contrary, the corrected BC (which is not consistent in the free 
transport regime) gives results that are far from the other schemes close to the boundary: 
it gives a numerical boundary layer at the left boundary that makes the solution much too 
large. 

5.3 Implicit diffusion 

Here, we illustrate the properties of the UGKS modified to recover an implicit scheme in 
the diffusion limit (see section I4~T]) . For this scheme, we take a time step defined by At = 
max(0.9eAx, c//Ax)-where cfl depends on the problem (between 0.9 and 0.1)-such that we 
get At = cflAx in the diffusion limit instead of At = cflAx 2 /2fi. 

In figure d we compare this modified UGKS (denoted by UGKS-ID in the following) 
to the non-modified UGKS (denoted by UGKS-ED) and to the reference solutions for two 
unsteady cases. The data are those of example 1 (with a kinetic regime) and 2 (with a 
diffusion regime) shown in section 15.11 For example 1, we are in a kinetic regime where 
our modification is useless. However, with the fine mesh, the two UGKS give results that 
are very close. With the coarse mesh, there are some differences, but still quite small. For 
example 2, we are in a diffusion regime, where the UGKS-ID allows to take a time step which 
does not respect the parabolic CFL At < Indeed, in the case of the fine mesh, while 
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the UGKS-ID requires a time step of 3.37 1CT 5 , the UGKS-ED requires a time step of only 
4.5 10~ 2 , which is 133 times larger. In addition, we observe that both schemes give results 
that are very close. Of course, the results are different for the small times (the schemes are 
first order in time, hence a larger time step induces a larger numerical error), but we see 
that for longer times, both schemes give very close results. As expected, if the time step of 
the UGKS-ID is decreased, the difference with the UKGS-ED reduces too. 

In figure we compare our schemes for long time cases. The data are those of examples 
3 (see section 15.11) . and 4 (see section 15.21) . We observe that for the coarse mesh, both 
schemes are very close. For example 3, the time step of UGKS-ID is 25 times as large as for 
UGKS-ED, while it is 16 times as large for example 4. For the fine mesh, UGKS-ID requires 
cfl = 0.1 (it is unstable for larger CFL), and the time step is here 100 times as large as for 
UGKS-ED in example 3, and 10 times as large as for example 4. For example 4, UGKS-ID is 
different from the reference solution, but this difference is as large as the difference observed 
for UGKS-ED. 

This study shows that the UGKS-ID can be efficiently used when the parabolic CFL 
condition is too much restrictive, with the same accuracy as the non- modified UGKS-ED. 

6 Conclusion 

In this paper, we have shown that the unified gas kinetic scheme, originally designed for 
rarefied gas dynamics problems, can be applied to other kinetic equations, like radiative 
transfer models. Moreover, the UGKS has been shown to be asymptotic preserving in the 
diffusion limit of such equations, as well as in the free transport regime. This scheme 
turns out to be an efficient multiscale method for kinetic problems, with a wide range of 
applications. In addition, we have shown that the UGKS can be simply modified to account 
for boundary layers and to obtain an implicit scheme in the diffusion limit. Finally, we have 
suggested an extension of the UGKS to general collision operator (in a non relaxation form) 
like that of neutron transport, that is still to be tested. This extension preserves the AP 
property of the method under some conditions. 

It would be interesting to rigorously prove the stability of the UGKS, in particular to 
determine an explicit CFL condition: we believe that the energy method we proposed for 
another AP scheme in [29J could be applied. For practical applications, we would like to 
extend this work to multidimensional problems. This should not be difficult, since it has 
already been done in rarefied gas dynamics by Xu and Huang in [TO] . 
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Figure 1: Kinetic regime (example 1): comparison between a reference solution and the 
UGKS (25 and 200 grid points). Results at times t = 0.1, 0.4, 1.0, 1.6 and 4 (e = 1). 
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Figure 2: Diffusion regime (example 2): comparison between diffusion solution and the 
UGKS (25 and 200 grid points). Results at times t = 0.01, 0.05, 0.15 and 2 (e = 10~ 8 ). 
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Figure 3: Intermediate regime with a variable scattering frequency and a source term (ex- 
ample 3): comparison between a reference solution and the UGKS (40 and 200 grid points). 
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Figure 4: Boundary layer problem (example 4): stabilized, corrected, and blended boundary 
conditions for the UGKS are compared to the explicit scheme (solid line with boundary 
layer) and the diffusion solution (solid straight line), e = fCT 2 . 
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Figure 5: Boundary layer problem (example 5): stabilized, corrected, and blended boundary 
conditions for the UGKS are compared to the diffusion solution, e = fCT 4 . 
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Figure 6: Free transport regime (example 6): stabilized, corrected, and blended boundary 
conditions for the UGKS are compared to the upwind explicit scheme, e = 1, a = 0. 
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Figure 7: Comparison of the UGKS with/without implicit diffusion with a reference solution 
for unsteady cases. Data of of examples 1 (kinetic regime, top) and 2 (diffusion regime, 
bottom). 
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Figure 8: Comparison of the UGKS with/without implicit diffusion with a reference solution 
for long time cases. Data of of examples 3 (Intermediate regime with a variable scattering 
frequency and a source term, top) and 4 (intermediate regime with non-isotropic boundary 
conditions, bottom). 
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